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ABSTRACT 

The dense molecular cloud cores that form stars, like other self-gravitating objects, 
undergo bulk oscillations. Just at the point of gravitational instability, their funda- 
mental oscillation mode has zero frequency. We study, using perturbation theory, the 
evolution of a spherical cloud that possesses such a frozen mode. We find that the cloud 
undergoes a prolonged epoch of subsonic, accelerating contraction. This slow contrac- 
tion occurs whether the cloud is initially inflated or compressed by the oscillation. The 
subsonic motion described here could underlie the spectral infall signature observed in 
many starless dense cores. 

Subject headings: ISM: clouds, kinematics and dynamics — stars: formation 



1. Introduction 

The central process in star formation is the gravitational collapse of a dense, molecular cloud 
core. Such objects are well studied empirically, through a variety of techniques. Less than half of 
observed dense cores already contain embedded stars, with the fraction increasing at higher mean 



column density and volume density (jBeichman et al 



1986 



JessoD Ward-ThomDSonll200m . The 



remai ning, so-called star l ess, cores exhibit a range of density contrasts and temperature gradients 



(e.g., iTafalla et al 



2004 



Crapsi et al.l 120071). Those of higher density do not differ in their gross 



properties from cores with stars (jLee &: Mverall999l ). Observations, then, have not yet provided a 
clearcut answer to the salient question: How does the collapse process actually begin? 

One important clue is the low level of internal motion observed within dense cores. This 
fact was esta blished early, from the relatively narrow linewidths of optically thin molecular tracers 
such as NH3 (jMvers &: Bensonlll983l ). The linewidths reveal gas motion that is lar gely subsonic, in 
contrast to the turbulent or wave- like character of the cores' external environments (Goodman et al 
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19981 ). Other studies bolster the view from spectroscopy. lAlves et al.l (j200ll ) used near- infrared 
extinction mapping to reconstruct the density structure of the starless core B68. After azimuthal 
averaging, the density profile is close to a theoret i cal one for a n isothermal sphere in balance 
between self-gravity and gas pressure (jEbertl Il955l : iBonnoii Il956l ). Thus, the object appears to 
be in dynamical equilibrium, as expected when internal motion is subsonic . Similar results hav e 
since been found for a number of other starless cores (jTafalla et al.l 12004 : iKandori et all l2005l ) . 
The inferred de nsity profiles are consistent w ith earlier ones reconstructed through submillimeter 
emission maps (IWard-Thompson et al.lll999l ). 



Real dense cores are not perfect spheres, in part because they d erive additiona l support from 
the anisotropic force associated with the interstellar magnetic field (jCrutcheii Il999l ) . In any case, 
the fact that clouds begin collapse from a near-equilibrium state is significant theoretically. What is 
the nature of this equilibrium? Clearly, the initial state cannot be dynamically stable with respect 
to small perturbations, as no collapse would ensue0 On the other hand, it is also unlikely to be 
dynamically unstable; it is difficult to imagine the prehistory of an object that could have veered 
away sharply from equilibrium at any time. The most natural initial state just prior to collapse 
for a star-forming dense core is one of marginal stability. While the object is in force balance, its 
fundamental mode of oscillation has zero frequency. Working again within the idealized framework 
of isothermal spheres, our paper addresses two specific questions. How does a marginally stable, 
isothermal cloud evolve with time? Can we account in a simple and compelling way for the tendency 
of such an object to collapse? 

Collapse calculations have a long history, and we are by no means the first t o recognize th e 
particular importance of t h e marginally s table isothe rmal sphere a s an in itial state. iHunteii (|l977l ). 
Foster fc Chevalien (|l993l ) , lOgino et al.l (|l999l ) , and lAikawa et all (|2005l ) all followed the collapse 
of such a configuration. However, these authors were concerned with developments at a relatively 
advanced stage. The subtle question of how the cloud first evolves away from equilibrium was 
effectively bypassed through assuming that the cloud was initially slightly overdense with respect 
to equilibrium, thereby guaranteeing eventual collapse. Our calculation focuses entirely on the 
transition issue. We show, using the tools of perturbation theory, that marginally stable clouds 
enter a protracted phase of slow, but accelerating inward contraction. We do not follow the evolution 
deep into collapse, a task that could be pursued using more standard methods. 

Our study of this early contraction phase holds more than pur ely theoretical interes t . A large 



fract i on of starless cores show convincing signs of inward rnotion (jWilliams et al 



199S 



Gregersen fc Evanal2000l : iLee et al 



2001 



199S 



Lee et al 



Schnee et al.ll2007l ). The spectroscopic signature of 



infall is an asymmetric emission line, often self-absorbed, that is skewed toward the blue. Inferred 
speeds are well below thermal, but greater than those associated with gravitational settling me- 



^ Starless cores of the l owest den sity contrast are indee d stable, and are confined by external pressure more than 
self-gravity. Andre et al. 1 2008h and Keto fc Casellil 1 2008h have suggested that many such objects are fated never to 
collapse, but will gradually disperse. 
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diated by ambipolar diffusion of the magnetic field (as calculated, e.g., by ICiolek &: Mouschovias 



19941 ). Our accelerating contraction is an attractive candidate to explain this intriguing, and per- 



haps pivotal, set of observations. 

In Section 2 below, we present our method of solution. After nondimensionalization of the 
dynamical equations, we introduce a perturbation expansion that allows us to separate out the 
first-order, oscillatory motion from the second-order displacement that progressively evolves with 
time. Sections 3 and 4 are devoted to both analytical derivations and numerical results. In the first, 
we display the cloud's normal modes of oscillation, including the critical one of zero frequency. In 
the second, we trace the bulk contraction of a cloud subject to the zero-frequency mode. Finally, 
Section 5 discusses the possible connection to the infall signature of starless cores, and indicates 
fruitful extensions of this work. 



2. Solution Strategy 
2.1. Physical Assumptions 



Our task of following early cloud evolution is facilitated by adopting a simplified physical 
picture. These simplifications are by now traditional, but it is important to revisit them as new 
observations arise and theoretical understanding grows. Thus, our employment of an isothermal 
equation of state, with a gas temperature that is constant in both space and time, is technically 
inconsistent with recent observational studies of starless cores. The inferred temperature of L1544, 
a starle ss core in Taurus t hat exhibits infall, decreases from 12 K at the outer edge to 6 K near the 
center (jCrapsi et al.l 120071 ). The central temperature is depressed by the partial shielding of dust 
grains from external starlight; these grains are thermally coupled to the gas through collisions. 
However, the region that is effectively shielded comprises little of the total mass, and collapse 
calcul ations that account for the finite temperature gradient show it to have a relatively minor 
effect (|Keto fc Field! bood v 



Another traditional assumption, and one ostensibly even more radical, is that our model cloud 
is spherical. It has long been known that the projected molecular-line emission contours of starless 
core s, like dense cores in general, are more accurately elliptical, with mean axis ratios of about 
2:1 (iLee fc MyersI Il999l ). The deprojected, or intrinsic, distribution of shapes must be inferred 
statistically from observations of a large sample. Under the condit ion that the under l ying structur e 



be axisymmetric, a prola ,te configuration best matches the data (iMyers et al 



also Jones et al 



2001 



1991 



Rvdenlll99ffl . 



Rela xing this restrictio n, I Goodwin et al.l ()2002l ) found that a triaxial configuration is preferred (see 



This result is puzzling. As we stated earlier, a departur e from sphe rical 
symmetry could partially be explained by additional magnetic support. However, iGallil (|2005l ) has 
demonstrated that there are no triaxial magnetostatic structures, at least for scale- free configura- 
tions. It may be, as Galli suggests, that the observed shapes reflect distortion created by internal 
oscillations, a phenomenon we ourselves shall invoke shortly. In any event, our adoption of a spher- 
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ical geometry, along with our neglect of magnetic support, can only be viewed as a convenient, first 
approximation. The essential finding of accelerating contraction from marginal stability should 
continue to hold within a more complete description of the equilibrium, magnetostatic state. 

Within the last decade, a number of researchers have questioned the very existence of equilib- 
rium structures. They have made the assertion, based on numerical sim ulations, that dense cores 



are tr ansient, fully dynamic entities. Such studies, recently reviewed bv iBallesteros-Paredes et al 



(|2007l ). treat an interior portion of the larger, parent cloud as a computational box filled with gas, 
with or without an embedded magnetic field. Turbulence is simulated by stirring the gas. This 
energy input is eventually offset by the unavoidable numerical dissipation present in the calcula- 
tion. Once such a steady state is attained, external stirring is turned off and self-gravity turned on. 
Sufficiently overdense structures collapse on themselves, in apparent imitation of the star formation 
processEI 

It is intriguing that these structures resemble, in their masses, sizes, an d geometric aspec t 
ratios, observed dense cores. Even triaxial structures have been obtained (e.g.. lOffner et al.ll2008l ). 
However, the key observational distinction between dense cores and their surroundings is their 
subsonic (or, more properly, sub-Alfvenic) internal velocity. This fact implies, as we have stressed, 
that the entities, unlike their numerical surrogates, are in force balance. Moreover, the large 
observed fraction of starless core s indicates a correspondingly protracted epoch in a typical dense 



core' s lifetime prior to collapse (iLee &: Myers 



)rresp 
I999I : 



Jessop fc Ward-Thompson 



200c 



Kirk et al 



2OO5I ). If these structures indeed condense from a turb ulent medium, their initial growth may be 
retarded by a relatively strong ambient magnetic field (jGalvan-Madrid et al.l 120071 ). 



2.2. Dynamical Equations and their Nondimensionalization 

We consider a spherical cloud if fixed mass M and isothermal sound speed a. The equation of 
state relating the internal pressure P and density p is 

P = pa"" . (1) 

The outer edge of the cloud is set by the condition that P falls to some external value Pe, also 
constant. In spherical symmetry, any matter outside the boundary exerts no gravitational force on 
the cloud. Hence we do not need to consider this medium in detail, beyond the fact that it exerts 
a fixed pressure. In particularly, we need not assume (unrealistically) that the external gas has 
a relatively high temper ature and low density, as is frequently done in collapse calculations (e.g.. 



Foster fc Chevalieijll993l ). We shall see that the cloud boundary, as just defined, shrinks during 



the evolution. This finding suggests that the distinction between "cloud" and "external medium" 



^In some simulations, t urbulence is continually driven, even after self-gravity is turned on (e.g., 
Vazquez-Semadeni et al. 20051 '). 
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is somewhat artificial, and that inward motion is spatial ly extended. Inde ed, the observed infall 
signature of starless cores is striking for its broad extent (jMyers et al.ll200Cl ). 



In Eulerian coordinates (r, i), the equation of mass continuity is 

dp _ 1 d{r'^ pu) 
dt dr 

where u{r, t) is the fluid velocity. This velocity satisfies momentum conservation: 

du du a? dp Gm 
dt dr p dr r 

Here, m is the mass enclosed within any radius: 



(2) 



m = Airr'^ pdr . (4) 
Jo 

As the total cloud mass is fixed, we shall find it more convenient to switch to Lagrangian 
coordinates {m,t). The mass continuity equation (2) becomes 

dp . 2 d f 2dr\ 

di = -''''p d^V di) ^ 

Here, we have replaced the velocity u by the time-derivative of the (now dependent) variable r: 

dr 

and have again suppressed all subscripts in the partial derivatives. The momentum equation (3) 
transforms to 

— = _47rr2a2-^ - ^ (7) 



Finally, equation (4) for m is replaced by 

dr 



dm Anr'^ p 



(8) 



Equations (5), (7), and (8) constitute our basic dynamical equations. We may further cast 

them into nondimensional form using the three constants G, a, and Pe- We define nondimensional 
versions of m, r, p, and t: 

m = J , (9) 

r = 2 ' (10) 

P - ^. (11) 

I . iil^. (12) 
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We further replace p by a new nondimensional variable (j), defined through 

~p = e-K (13) 
After dropping the tilde, our dynamical equations become 



d(t> 
'dt 


dm \ dt J ^ 


(14) 


w 


dm ' 


(15) 


dr 
dm 




(16) 


We will also be using equation (6) relating 


u and r. If we set u = 


u/a, then this equation remains 



the same nondimensionally. 



2.3. Perturbation Expansion 

Equations (14)-(16) are to be solved subject to the inner boundary condition r{0,t) = and 
the outer one (f){M,t) = 0. The latter is just the requirement, expressed in nondimensional lan- 
guage, that the gas pressure at the edge equal the fixed value Pg. We must also specify the cloud's 
initial configuration. This state is itself slightly perturbed from true equilibrium, under the influ- 
ence of the oscillation mode. We are thus motivated to introduce a perturbation expansion of the 
dependent variables about the equilibrium state. We write: 

<p{m,t) = (/)o(m) + e(/.i(m,t) + e2^2(m.,t) + ... , (17) 
r{m,t) = ro(m) + eri(m,t) + e^r2{m,t) + ... . (18) 

Here, e is a small, free parameter. The physical significance of e is that it will set the amplitude of 
the oscillation mode present at the outset (see §4 below). Quantities with the subscript refer to 
the time- independent, equilibrium state. 

We are now in a position to separate out the hierarchy of cloud motions. We substitute our 
series expansions into equations (14) - (16) and equate like powers of e. At the lowest (zeroth) 
order in e, both sides in the equation of mass continuity vanish. However, equations (15) and (16) 
yield, to the same order, the structural equations of the underlying, equilibrium state: 

(19) 
(20) 

When we equate terms proportional to e, all three equations give a non-trivial result: 

- frf-^) , (21) 



d(f)Q 


me"^' 


(III/ 




dro 




(III/ 


r2 



dt 5m V " dt 
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- ^ dm^^'^o'^^ ^o^ij+^oe + ^3 , (22) 

P - ^(^1-^) . (23) 
am ^0 \ ''o / 

These zeroth- and first-order equations are sufficient to describe small oscillations about a 
general equilibrium state. However, we arc interested in that unique equilibrium for which the 
fundamental oscillation frequency vanishes. The perturbed cloud will not be perfectly static, but 
will evolve slowly, through second-order effects. Thus, we also need to equate terms in the dynamical 
equations proportional to e^. We find 

, d(hi 862 rh d f ^ dri 9 dr2\ 



e-<^ 






dm \ 


Ai 


+ A2 : 




fSrf 







(25) 



p = !^ _ + ^ _ + I ) . (26) 

dm ri \ ri ro ro 2 ' 



The two acceleration terms Ai and A2 in equation (25) are 

A, = e-^°^(^-r2</.2 + r?-2ron0i+2ror2) 
— (2ron-ro0i)+roe — 



(27) 

A. d(h\ , o \ d(l)2 2 m r-i Smr? , ^ 

= «"'°^(2-ori-ro20i)+ro2e-^°^ + — ^. (28) 



3. Cloud Oscillations 
3.1. Equilibria 

We consider first equations (19) and (20), that describe the equilibrium state. These equations 
are subject to the boundary conditions ro(0) = and (f)o{M) = 0, to be applied at the center 
and outer edge, respectively. More accurately, the second condition defines the boundary of the 
equilibrium configuration. As is well known, there exists a one-parameter family of such structures, 
each characterized by its own center-to-edge density contrast. We denote this contrast as pc, since 
it is identical to the nondimensional central density in our formulation. 

Because of the inner boundary condition, the righthand sides of both equations are singular at 
the origin. To start our numerical integration, we therefore expand both ro(m) and 4'o{m) in the 
appropriate power series: 

ro(m) = aom^/^ + aim + 02 m^/^ + ... , (29) 
0o(to) = -Inpc + bom^/^ + bim'^/^ + b2m^ + ... . (30) 
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The first two coefficients are 



do 



bo 



1/3 



o 4 ' 
2pc 



(31) 
(32) 



and higher ones are determined recursively. 



Based on our numerical integration, the dashed curves in Figure 1 are density profiles for states 
with pc = 5 and pc = 30. Note that we have plotted po = exp(— </)o) as a function of tq, in the 
conventional manner. The middle, solid curve corresponds to pc = 14.04, for which M = 4.19. 
This is the famous Bonnor-Ebert state. As we shall verify shortly, it is the starting configuration 
of interest for the present problem. The other two cu rves represent outer limits to the empirical 
fitting of density profiles for most starless dense cores (jKandori et alll2005l ). 



3.2. Normal Modes 

Equations (21)-(23) describe, to linear order, internal motion of the equilibrium cloud. All 
such motion can be decomposed into a series of normal modes. The fundamental, or breathing, 
mode is generally expected to have the largest amplitud e, but higher harmonics may also be present 



(as well as non-radial oscillations; see lKeto et al.ll2006l ) 



To obtain the full set of spherical normal modes, we first solve equation (23) for the density 
perturbation: 

0i = ro^e-^°|l + ^. (33) 
om ro 

We substitute this expression into the momentum equation (22) and obtain a partial differential 
equation for ri(m,t): 

-oe-'^° ^ + 2e-^° (2ro - m) + ( - -A n . (34) 



dt"^ dm? dm \ r"^ ^ 

The coefficients in this equation are all functions of the chosen equilibrium state. 

Before proceeding with the solution of equation (34), we note that the additional first-order 
equation (21) has not been used in the derivation. In fact, this relation is redundant, and already 
contained implicitly in equation (33). To see this, take the time derivative of the latter: 

2 dri , , 

+ - ^ , 35a 
ro at 



dcPi 




dt 


" dt dm 








^ dt dm 
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which is just equation (21). Here we have used the time-independence of all equilibrium quantities, 
and have also employed equation (20). 

Returning to equation (34), we seek oscillatory solutions for ri(m,t). Thus, we set 

ri(m,t) = A Ciim) e'"^^ . (36) 

Since equation (34) is linear in ri, the coefficient A is an arbitrary constant. After substitution, we 
find an ordinary differential equation for ^{m): 



(37) 



The boundary conditions are now ^i(O) = and <j)i{M) = 0. Using equation (33), the latter may 
be transformed into a condition on ^i(m): 

d£,i 2^16^° 

— = o — ■ at m = M (38) 

dm Tq 

We note that equation (37) has a regular singular point at the origin. Once more, we start the 
integration through a power series development: 

^i(m) = com^/^ + cim + C2m^/^ + ... . (39) 

The constant cq may be chosen arbitrarily. The next coefficient is 

ci = . (40) 

and successive ones may be similarly found. 

For any central density pc of the equilibrium cloud, there exists a sequence of w^-values such 
that ^i(m) meets the two boundary conditions. When pc = 14.04, the lowest value of u'^ is zero. 
The corresponding ^i(m) is the oscillation mode of interest. If we denote the fundamental, first, 

and second harmonics as ^iq, ^h, and ^12, then the corresponding a;^-values are 0, 8.37, and 24.3. 
Figure 2 plots these three normal modes, again using ro as the independent variable. In all cases, 
we have set the coefficient cq in equation (39) to unity. 

Note finally that we may use the equilibrium relations, equations (19) and (20), to rewrite 
equation (37) as 

This equation is of the Sturm-Liouville type. As expected physically, all of its eigenvalues a;^ are 
real. Distinct eigenmodes are orthogonal: 



M 

iip^iqdm = . for p ^ q (42) 
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Since only the zero-frequency oscillation will be of interest in the following discussion, we will 

soon revert to the simpler notation ^i(m) for the fundamental mode, and assume the standard 
normalization cq = 1 for this function. After also setting the coefficient A in equation (36) to 
unity, we will thus be making the identification ri{m,t) = ^i(m). 



4. Cloud Contraction 

4.1. Fundamental Equation 

Once we have selected the oscillatory mode of the equilibrium state, the second-order equa- 
tions (24)-(28) describe additional motion. This motion is, of course, critical when the underlying 
oscillation has zero frequency, but we shall first keep the discussion more general. As in our deriva- 
tion of the normal modes, it is convenient to solve for the density perturbation (now ^2) from the 
mass-radius relation, equation (26): 

dr2 2 r2 rf Tq e~^ f dri 



^0 ® + ^2 2 \ dm / ' ^ 

where we have used equation (33) for (f)i. Wc substitute both this latter relation and equation (43) 
for 02 into the righthand side of the momentum equation (25). We thus derive our fundamental 
partial differential equation for the displacement r2(m,i): 

= rie-^'- is + 2e^*. (2.„ - „) |5 + (ij - |) + F , (44) 



dt'^ ° dm? dm 



where 



3 m\ 2 -w,- /2m 



rl 



-2</,o i2T.-m) 



F^[-,-^]ri + ^-'^[—-'^]^ + 4e-^''°{'^r,-m)[^\ . (45) 



Wc again remark that it has been unnecessary to invoke the continuity equation (24). The 
explanation, as before, is that this equation yields no new information. To see this, take the time 
derivative of equation (43) and multiply through by e*^" : 

e^o ^ ^ ^2 ^ 2e^ dr2 _ 2ri c'^° djn _ ^4g-^o ^ ^'^i _ (4g) 

dt ° dtdm ro dt dt ^ dm dtdm 

After using equation (33) for (f>i and equation (35c) for dcpi/dt, we also have 

2 dri dri ^ 4 dri d'^ri ^ ^ dro dri ^ ^ ^^^1 (47b) 
dm dt ^ dm dt dm dm dt dt dm ' 
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where we have again utiUzed equation (20). Combining equations (46) and (47b) then gives 



which is equivalent to equation (24). 

Except for the inhomogcncous term F, equation (44) for r2{m,t) is identical to equation (34) 
for ri{m,t). It is precisely because of this extra term that the second-order motion cannot be 
described as a normal mode of oscillation. Once we single out the zero-frequency, fundamental 
mode for the first-order motion, F can be written purely as a function of m: 



(49) 



The boundary conditions on r2{m,t) are the usual ones: r2{0,t) = and <^2{M,t) = 0. Prom 
equation (43), the outer boundary condition is more usefully recast as 



dr2 
dm 





/3r? 


2r2 






ro 


Q<P0 


(Hi 


2r2 









(50a) 

at m = M (50b) 



4.2. Method of Solution 

Solving the partial differential equation (44) requires that we first specify the full initial state 
of the cloud. Assuming the fundamental mode has zero frequency, i.e., that ri = we still need 
to set the functional form of r2(m,0). As our fiducial initial state, we demand that the cloud's 
density profile be just that resulting from the normal mode acting on the equilibrium state. That 
is, we set <^2{'m,0) = 0. With this condition, equation (17) then gives the physical interpretation 
of the parameter e as the nondimensional amplitude of the initial (^perturbation. 

The vanishing of (l)2{m,0) does not mean that r2(m, 0) = 0. Instead, the initial r2-profile 
follows by setting (f)2 = in equation (43) and specializing to ri = ^i: 

We solve this ordinary differential equation for r2(m, 0). After again noting the regular singular 
point at the origin, we begin the numerical integration through a power-law expansion: 



r2 = do m^' + dim + d2 m^'-^ + 



(52) 
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The expansion coefficients are readily found: 

do = . etc. (53) 
zao 

Tfie top curve of Figure 3 displays tfie calculated r2(m, 0), again as a function of tq. 

We also need to specify the initial velocity dr2/dt (m, 0). An interesting, and physically relevant 
situation is when the cloud is perfectly static, i.e., when 

<-) 

More realistically, the cloud has slowly evolved from some earlier state, perhaps through the effect 
of ambipolar diffusion. In the absence of a more complete evolutionary picture, we shall adopt the 
simplest, zero- velocity, initial condition. 

Equation (44) may be solved numerically through the method of characteristics. For this 
purpose, it is convenient to adopt ro as the independent, spatial variable. Using the connection 
between ro and m in equation (20), the fundamental equation (44) becomes 

+ {-- -2)1^ +{^- -2)^2 + (55) 



dt"^ dr^ \ro r^J dro 
where F is now written as 



2 3mV2 _^ /2m 2 \ d^f ^ / 2 m\/d6^^ 

The reason for changing Lagrangian variables from m to ro is that the characteristics of equa- 
tion (55) are simply 

^ = ±1. (57) 
dt ^ ' 

In light of our nondimensionalization, we see that small disturbances in the cloud travel at the 

isothermal sound speed, an intuitively appealing result. We solve equation (55) by propagating 

the partial derivatives dr2/dro and dr2/dt along the (ro,t) grid, as detailed in the Appendix. The 

starting value of dr2/dt is everywhere zero, according to equation (54). As we advance in time, we 

also incorporate both the inner boundary condition ro = and the outer one of equation (50b). 

The latter is now more conveniently written as 



dr2 _ 3Cf 2r2 
dro ro ro 



2 . at m = M (58) 



4.3. Numerical Results 



Integration of the partial differential equation (55) for r2(ro, t) shows that this function globally 
decreases for t > 0. That is, all mass shells contract, although the cloud was neither expanding nor 
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contracting in its initial state. If this state were slightly compressed from equilibrium, the parameter 
e would be negative. If slightly expanded, e would be positive. In either case, the second-order 
effects governing subsequent evolution are proportional to e^. Thus, the cloud inevitably contracts. 

This contraction begins slowly, but accelerates with time. To illustrate the effect, Figure 3 
displays r2(ro,t) for t = 0, 1, and 2. We see how any initial, slight inflation (assuming e > 0) is 
everywhere reversed by t = 1. By t = 2, mass shells roughly midway from the center are clearly 
contracting the most rapidly. 

The time span shown in the figure exceeds that associated with free-fall collapse. The latter 
is conventionally taken to be 

where pc is the cloud's central density. (Here, we are temporarily reverting to dimensional variables.) 
Our nondimensional time t is given by equation (12), so that 



- = -J^t. (60) 

Here, pe = Pe/o?- For the Bonnor-Ebert initial state {pc/pe = 14.04), the numerical coefficient 
in equation (60) is 1.95. Thus, the last profile shown in Figure 3 corresponds to 2 x 1.95 = 3.90 
free-fall times. 

We are primarily interested, not in the displacement of mass shells, but in the contraction 
speed. From equation (18), the nondimensional velocity (i.e., its value relative to the sound speed 
a) is given by e^dr2/dt. The ac tual velocit y of ea ch mass shell thus depends on the value of e. 



which we do not know a priori. iKeto et al.l (|2006l ) have recently matched molecular-line profiles 



from the starless core B68 by assuming the object is undergoing a nonradial oscillation mode, 
with a dimensionless amplitude of 0.25. Using this figure as rough guide, we provisionally choose 
e = -1-0.2 and display the resulting velocity profiles in Figure 4. Here, the times are identical to 
those in Figure 3. However, the independent spatial variable is now the full radius r, as given by 
equation (18). 

We see that the velocity of all mass shells is negative, with a magnitude that increases with 
time. Thus, contraction is indeed accelerating. By the last time shown, the largest contraction 
speed within the cloud is about 0.2 times the sound speed. Interestingly, the cloud edge also has 
negative velocity, suggesting that contraction is spreading into the exterior region. 

Closer to the cloud center, it appears from Figure 4 that the velocity at any time increa ses 



linearly with radius. This feature of the evolution was found by iFoster Chevalieii (|l993l ) in 
their numerical simulation. We may derive the result analytically by examining the fundamental 
equation (55) in this limit. Consider first the inhomogeneous term F. Now the series expansion 
for ^i{m) in equation (39) tells us that = eo^o to lowest order, where cq is a constant. From 
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equation (56), F has three terms that diverge as tq approaches zero. Their sum is 

Mi _ 1 M + 1 ^ 2e| _ 4e| ^ 2e| ^ 

rl rl dro ro \dro J ro ro ro 

which vanishes. However, the rest of equation (55) has coefficients that still do diverge. Close to 
the origin, the equation reduces to 

d^r2 2 dr2 2r2 , . 

= + — a T ■ ro^O (62) 

dr^ ro dro 

Here we have used the fact that, since r2 = at the origin for all times, we may neglect the second 
time derivative on the left side of equation (55). Since equation (62) itself holds for all times, we 
may take its time derivative and obtain the analogous relation for the velocity v = e^dr2/dt: 

d\ 2 dv 2v 
= + J . ro^O 63 

The non-divergent solution to this equation is that v is indeed proportional to ro. 

Figure 5 shows the evolution of the cloud's density profile p(r, t) over the same time interval 
as in Figures 3 and 4. Since e was assumed to be positive, the cloud begins in a slightly inflated 
state, with a central density of only 5.15. However, subsequent contraction drives up the density, 
which reaches a central value of 15.5 by i = 2. Since all internal velocities are still subsonic at this 
time, the density profile is consistent with a cloud that, to first order, is in hydrostatic balance. 
Force balance will, of course, be strongly violated in the future, as the cloud enters a state of true 
collapse. 

The perturbative nature of our calculation limits us to describing the initial transition phase. 
Our calculation is only valid to the point where second-order terms in the expansions become com- 
parable to their first-order counterparts. For e = -1-0.2, the maximum absolute value of r2(m, t) 
in equation (18) becomes equal to the maximum value of e ri(m, t)att = 1.9. Thus, the calculation 
is still marginally valid for the t = 2 profiles shown in Figures 3, 4, and 5, but not beyond those. 



4.4. Alternative Initial Condition 

It is important to verify that the evolutionary results shown thus far are not sensitive to the 
detailed initial state chosen. Recall that our fiducial state represented a pure first-order pertur- 
bation of the density, in that we set 4'2 = ^ for all mass shells. This choice was convenient, but 
arbitrary. We could have chosen any initial density profile consistent with the boundary condition 

MM,t) = 0. 

Another starting configuration is obtained by setting 

2m) ■ ^^^^ 
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With this 02-pi'ofile, the outer boundary condition is satisfied, but the central (j)2 is now unity. 
Since p = exp (—(/>), the initial cloud is less dense than before. If we insert equation (61) into 
equation (43), we may again integrate an ordinary differential equation for r2(m,0). As shown in 
the left panel of Figure 6, the r2-profile is thus inflated relative to the previous initial state, with a 
maximum fractional difference of 43 percent. Nevertheless, subsequent contraction brings the cloud 
to a very similar configuration. The right panel of Figure 6 displays the new r2-profile at t = 2. 
It differs from the one obtained using the original starting state by at most 8 percent. The profile 
of velocity, dr^jdt^ is also very close to the previous one, with the maximum speeds differing by 
less than 1 percent at t = 2. 



Discussion 



In the picture introduced here, star-forming dense cores undergo a prolonged phase of con- 
traction before their ultimate collapse. The contraction is slow because the cloud's self-gravity is 
still nearly counteracted by the outward pressure gradient. The slight imbalance of these forces 
creates subsonic, inward motion that gradually accelerates. Previous researchers performing direct 
simulations have also documented slow motion prior to runaway collapse. However, the character- 
istics of this phase, and indeed whether it led to collapse or rebound, depended on the unavoidable 
artifices of a direct simulation: the imposed overdensity in the cloud , the numerical accuracy of 



the initial state, and the precise treatment of th e central few zones (jHunteri 119771 : iBoss &: Black 



19821 : [Foster &: Chevalieij Il993l : lOgino et al.l |l999j) . We have demonstrated physically how an ini- 
tially static, marginally stable cloud inevitably contracts. This contraction is essentially the inward 
phase of a slow oscillation that smoothly leads to free-fall collapse. 

Motivated by the observed infall signature of many starless dense cores, others have offered 
different pictures. iKeto &: Field! (120051 ) hypothesized that some cores are born in a gravitationally 
unstable state. They followed the collapse of such an object numerically, accounting for internal 
cooling by molecular lines and thermal dust emission 



This work extends that of Zhou et al 



(|l99d l , who modeled the (starred) core B335 as undergoing collapse from a singular isot hermal 
sphereo As already noted, the early history of unstable objects is proble matic; Keto &: Field (2005) 
speculate that the cloud fragmented from a larger, turbulent velocity field. iMvers &: LazarianI (|l998l ) 
attributed this localized fragmentation to enhanced dissipation via ion-neutral friction, leading to 
a pressure-driven, inward flow. 

Our general concern about using unstable or actively collapsing states to match observations is 
their brevity. To illustrate the point more quantitatively, suppose the infall signatures reflect first- 



The si nsular isothermal sph ere is unstable not only to the fundamental oscillation mode, but to all higher 



harmonics (Stabler &: Palla 



20041 . Chapter 9). The collapse of an u nbounded, sin g ular s phere was calculated in a 



self-similar fashion by Shul 1 197^ . In light of the infall observat 
include a finite, inward velocity at large radii 



Fatuzzo et all l|2004l ) generalized this model to 
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order motion, i.e., that the cloud's fundamental eigenfrequency is nonzero. The nondimensional 
perturbation ri(m,t) would then be expressed as 

n = -^1 e+l'^°'* . (65) 

The function ^i(?ti) obeys equation (37), but with replaced by — jwoP- If we again normalize the 
coefficient Co in equation (39) to unity, then ^i(m) resembles the curve ^lo in Figure 2. Notice the 
overall minus sign in equation (65) that ensures contraction prior to free-fall collapse. Over what 
time does this transition occur? 

The physical velocity, in units of the sound speed, is edri/dt = e |a;o| ri. If e is still about 0.2, 
and if we are to match velocities of 0.2 times the sound speed at the present epoch [t = 0), then 
the growth rate |wo|~^ must be about unity. The dimensional time t for the velocity to increase by 
a factor e is, from equations (9) and (12), 

GM , , 

t = . 66 

The numerical solution to the altered equation (37) tells that the nondimensional cloud mass M 
corresponding to |a;o| = 1 is 4.02. The center-to-edge density contrast of this object is 30, a figure 
that is margin ally consistent with obs ervations (recall Fig. 1). Returning to L1544, its gas temper- 



ature is 10 K ([Benson &: Mverslll989l ). while the most recent mass estimate is 2 Mq (jOhashi et al 
19991 ). Equation (66) then gives an e- folding time of 2 x 10^ yr. The statistical prevalence of 
starless dense cores makes it unlikely that they are evolving over such a brief, dynamical interval. 
Note finally that if the initial perturbation amplitude e were smaller than 0.2, then |u;o| would be 
correspondingly larger, driving down the evolutionary time even more. 

We emphasize again the general nature of our theoretical finding. A more realistic model 
for a starless core should certainly include the anisotropic supporting force from the interstellar 
magnetic field. Such magnetostatic configurations are themsel ves subject to global o scillations. The 



lowest eigenfrequency vanishes in the marginally stable state (jTomisaka et al.lll988l ). These states, 
analogues of the Bonnor-Ebert configuration used here, will also undergo slow contraction prior 
to collapse, even under the approximation of flux freezing. They are close to bei ng magnetically 



supe rcritical, so that ambipolar diffusion may enhance the contraction process (jCiolek &: Basu 



200 ll ). In any event, it will be interesting to track the evolution through the transonic phase into 



full collapse, both in our spherical model and its magnetized generalization. 

Returning to the observational motivation of this study, it will also be instructive to calculate, 
within our spherical model, the predicted profiles for molecular emission lines of varying optical 
depth. For our representative case of e = -1-0.2, we find a maximum infall speed of 0.2 times the 
sound speed after about 4 free-fall times. This speed is comparabl e to the 0.08 km s~^ inferred 



for the well-studied starless core L1544 through N2H"'" observations (jWilliams et al.lll999l ) We are 
encouraged by this finding, but stress the need for more comprehensive and detailed comparisons. 
Note especially that the amplitude e is not readily observable with any precision; nor is the evolu- 
tionary time t. And yet contraction models of lower e and larger t broadly mimic, in their velocity 
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profiles, those with higher e and shorter t. Hopefully, calculated line profiles will differ enough to 
resolve this ambiguity. If some of the profiles successfully match observations, we will not only have 
gained new insight into the mechanism of dense core contraction, but also a new measure for their 
longevity prior to collapse. 

We are grateful for stimulating conversations with Phil Chang and Steve Shore during the 
inception of this project. We also thank the referee, Tom Hartquist, for comments that improved 
the original manuscript. S. S. was partially supported by NSF Grant AST-0639743. 
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A. Implementing the Method of Characteristics 



Equation (55) is a linear, inhomogeneous parti al differenti al equation in the independent vari- 
ables ro and t. Following standard procedure (e.g., lAmeslll992l ). we suppose that along some curve 
in the (ro, t) plane, we know r2 and its first partial derivatives dr2/drQ and dr2/dt. Then the three 
second partial derivatives are related through equation (55), which we rewrite as 



52 r 2 92 



r2 



Here, 



9 



9^-+ ^ n]r2+ F 



2 m 



2 m\ dr2 
Jo r^J dro y '0 'o- 
and F is given by equation (56). Along our curve, we also have the differential relations 



92 r2 
dro dt 
92 r 2 



At 



At + 



d'^r2 
dr^ 
d^r2 
dro dt 



Aro 



Aro 



4'7r 

\dro 

a(^ 

V dt 



(Al) 
(A2) 

(A3) 
(A4) 



Equations (Al), (A3), and (A4) constitute three algebraic relations for the second partial 
derivatives. We may recast the system in matrix form 



where 



X 



y 



and 



M 



Mx = y , 

/ d'^r2/de 

d'^r2/dro dt 
\ d'^r2/drl 

( ' \ 

A{dr2/dro) 

\ A{dr2/dt) I 

1 -1 
At Aro 
At Aro 



(A5) 



(A6) 



(A7) 



(AS) 



The vector x is uniquely determined unless det M = 0. Thus, discontinuities propagate along 
characteristics given by 

Aro 



At 



±1 



(A9) 



as in equation (57). 
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We may also replace any column in M by the vector y. For example, consider the matrix M' 
given by 

/I . -1 \ 

M' = \ A{dr2/dro) Avq ■ (AlO) 
\ At A{dr2/dt) I 

There is no solution at all for the second derivatives unless detM' = 0. This condition gives us 
the further differential relation 



A(l^l-Ar^)=sA.o, (All) 



9ro / V 



along the + characteristic, and 



along the - characteristic. 



Figure 7 is a schematic spacetime diagram that illustrates the practical procedure. We set 
up uniform grids along the tq- and t-axes. At t = 0, we have dr^jQt = 0, while the solution 
to equation (51) gives the initial values of 5r2/5ro for the fiducial initial state. (An analogous 
equation is solved for the alternative state; see §4.4.) Starting from any two adjacent points on the 

ro-axis, such as A and B, we use equations (All) and (A12) to solve simultaneously for the two 
first partial derivatives at point D, where the two characteristics intersect. Similarly, we find the 
two derivatives at E starting from the pair B and C. The derivatives at D and E then yield those 
at point G, and so on. 

We also need to propagate the information contained in the central and surface boundary 
conditions. Since r2(0, i) = 0, it is also true that dr2/dt = at points such as F in the figure. 
This condition, along with equation (A12) for the - characteristic joined to point D, allows us to 
solve for dr2/drQ at F. Since both partial derivatives are now known at F and G, this information 
can then be used to find the derivatives at points further advanced in time. Near the cloud's 
outer edge, we establish the two partial derivatives at point L in the usual manner. Knowing this 
information at L, we find one relation between the two derivatives at the boundary point N by 
using equation (All) for the + characteristic. The outer boundary condition, equation (58) then 
supplies a second relation between the derivatives. To find r2 itself at N, we further use 
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Fig. 1. — Density profiles of equilibrium, isothermal clouds. Both the density and radius are in the 
nondimensional units defined in the text. The solid curve is the critical Bonnor-Ebert state, while 
the two dashed curves represent approximate outer bounds obtained by empirical fitting to starless 
dense cores. 




Fig. 2. — Normal modes of the Bonnor-Ebcrt isothermal sphere. Shown is the first-order displace- 
ment ^1 as a function of radius, for the primary (^lo), first harmonic (^ii), and second harmonic 
(^12). All curves have been normalized to have the same initial slope. 




Fig. 3. — Second-order displacement r2, shown as a function of the Lagrangian radius tq. Prom 
top to bottom, the three profiles correspond to t = 0, 1, and 2, respectively 




Fig. 4. — Nondimensional velocity profiles as a function of the full radius r. From top to bottom, 
the corresponding times are t = 0, 1, and 2. The initial state was perturbed with a dimensionless 
amplitude e = +0.2. 
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Fig. 5. — The evolving density profiles for a cloud perturbed initially with e 
density monotonically increases for the three times shown: t = 0, 1 and 2. 



= +0.2. The central 
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Fig. 6. — Evolution of the displacement r2 with two initial conditions. The left panel shows that, 
at t = 0, the r2-profile for the altered initial state (solid curve) differs substantially from that 
associated with the fiducial initial cloud that had a purely first-order density perturbation (dashed 
curve). Nevertheless, as the right panel shows, the profile is close to the previous one at t = 2. 
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Fig. 7. — Schematic spacctimc diagram illustrating the solution procedure for the fundamental 
partial differential equation (55). The straight, diagonal line segment from A to D lies along a 
+ characteristic, while the line segment from B to D lies along a - characteristic. 



